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ABSTRACT 

We construct a realistic model for the dark halo substructures, and derive analytically 
their spatial, mass, and velocity distributions. Basically, our model is a modification 
of the Press-Schcchter theory to account for the dominant dynamical processes that 
mark the evolution of dark halo substructures such as the tidal stripping and dynamical 
friction. Our analytic model successfully reproduces all the well known behaviors of 
the substructure distributions that have been found in recent numerical simulations: 
the weak dependence of the mass distributions on the host halo mass; the anti-bias 
of the spatial distribution relative to the dark matter particle components; the nearly 
power-law shapes of the mass and velocity distributions. We compare our analytic 
results with recent high-resolution N-body simulation data and find that they are in 
excellent agreement with each other. 
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1 INTRODUCTION 

The substructures of dark matter halos have attracted many 
attentions, especially because of a possible problem they 
poses in the currently popular cold dark matter (CDM) 
model. Numerical simulations have confirmed that the 
CDM model predicts roughly 10 % of mass in a dark halo 

is bou nd to subs tr uctur es llTormen. Diaferio fc Svep 

19981: iKlvpinetalJ Il999l: lOkamoto fc Habd Il999t 



Moore et alj Il999t iGhigna et alj l2000l: ISprineel et alJ 
20011: IZe ntner fc Bullock' '2003t iDe Lucia et al.l 120041: 

Kravtsov. Gnedin fc Klypin 2 00 . Although the current 
standard CDM paradigm has been very successful in 
matching the observational data on large scale, it has been 
noted that the CDM model overpredicts the abundance 
of substructures compared with that of observed galactic 
satellites (e.p;.. iKauffmann. White fc Guiderdonil Il993l : 
iMoore et all 119991) . This CDM problem on sub-galactic 
scale is regarded as one of the most fundamental issues that 
has to be addressed. 

It has been proposed that this CDM problem 
on sub-galactic scale could be resolved by changing 
radically the nature of dark matter. The propos- 
als include the sel f-interacting dark matter model 
JSpergel fc Steinhardtl .2000). the warm dark mat- 
ter model ilcolfa~Avila-Reese fc Valenzuelal I2OO0I : 
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iBode. Ostriker fc Turo3 1200 ll). the annihila, ting dark 
matter model |Kaphnghat. Knox fc Turneill2000D . and non- 
thermally produced dark matter (lLir^^ni20Mh . Another 
possibility is to introduce new inflationary models that can 
produce density fluctuations with small-scale power cut-off 
such as an inflation model w ith broken scale-invariance 
iKamionkowski_fcLiddlj 120001) and a double hybrid in- 
flation ifYbkovam?l200g) . However, these rather radical 
models are gradually in disfavor wi th recent ob servations 
and numerical simu lations (e.g.. iHennawi fc Qstrikeni20()3 : 
lYoshida et 'ai]l200^^ . 

The problem may be resolved also by taking account 
of astrophysical p rocesses such as photoionizing background 
JSomerviUd 120021) and inefficient star formation in small 
mass halos JStoehr et alJl2002D . These ideas claim that the 
observed number of satellite galaxies is small because only 
very massive substructures contain stars and most sub- 
structures are dark. To detect these small dark substruc- 
tures, the gravitational lensing can be a very powerful tool. 
For instance, we can constrain the amount of substruc- 
tures through anomalous fiux ratios [Ma o fc Schneider 199 ^: 
'Met calf fc Madaul l200lt IChibal bool balal fc Kochanekl 
I2OO2I) . s pectrosc'opy jMoustakas fc MetcalJ l2003l) . or mon- 
itoring llYonehara. UmernuK^^usn200a) . Although the 
gravitational lensing allows us to probe the distribution 
directly, a caveat is the result is dependent on the spa- 
tial d i stribution of substru ctures llChen. Kravtsov fc Keetoiil 
120031: lEvans fc Wittlf2003l) . However, these ideas based on 
the dark substructures may not be consistent with obser- 
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vations, either: Recent high-resolution numerical simula- 
tions have found that the massive substructures tend to 
place in the outer part of host halos, which is not the case 
for the satellites of our M ilky way dPe Lucia et alJ 120041 : 
iTavlor. Babul fc Silklliool . 

Despite the importance of substructures to understand 
the structure formation in the universe, most work had to 
resort to numerical approaches in studying mass and spa- 
tial distributions of substructures. Compared with the nu- 
merical approach, the analytic approach has an advantage 
that it allows us to compute distributions in a wide range 
of parameters, e.g., mass of host halos and power spectrum. 
This is essential in cosmological applications of substruc- 
tures, e .g., in computing power spectrum from halo ap - 
proach jSheth fc .Tain l l200,4 IPolnev. .Tain fc Tak adal l2004ll . 
In contrast, in N-body simulations one needs high force and 
mass resolutions to overcome "the overmerging" problem 
ijKlvpin et alJll99fll ). which makes the large number of reli- 
able calculations difficult. 

There are several attempts to derive the mass func- 
tion of sub structures by us ing analytic methods (e.g., 
iFuiita et al1 l2002: Shcthl l2003l') . These previous approaches, 
however, were all based on the oversimplified assump- 
tions such as the mass conservation and the uniform spa- 
tial distribution of dark halo substructures. They consid- 
ered only the gravitational merging of substructures, with- 
out taking account the effects of other important dynam- 
ical processes that drive the substructure evolution, lead- 
ing to their final distributions. Among them, most impor- 
tant are t he mass-loss caused by the global tides of host 
halos (e.g. lOkamoto fc Habelll999^ . and the orbital decay 
driven by the dynam ical frictions fe.g.. Ilbrmen et aLil99a: 



gynam ii 
]|2H). 



IKravtsov et al 

Very recently. iLed ll2004f) . for the first time, developed 
an analytical formalism for the global substructure mass 
function by incorporating the effect of tidally driven mass- 
loss effect and the non-uniform spatial distribution of sub- 
structures in dark halos. Yet Lee's approach was still a 
simplification of reality, using a crude tidal-limit approxi- 
mation in estimating the tidal mass-loss and also ignoring 
the effect of dynamical friction which may be more impor- 
tant in estimating the abundance of massive substructures 
IHavashi ct al. 2003). 

In this paper, we construct a more realistic model for 
the evolution of dark halo substructures than previous ap- 
proaches by taking account of not only the effect of global 
tides but also the effect of dynamical friction, and derive 
analytically the spatial, mass, and velocity distributions of 
substructures in dark halos. We also compare our analytic 
results with very recent high-reso l ution numerical simula- 
tions presented bv lDe Lucia et all i2004f) . 

The plan of our paper is as follows. In 321 we study the 
effects of global tides and dynamical frictions on the evolu- 
tion of dark halo substructures, and provide approximation 
formula to quantify them. Section |3] is devoted to present 
an analytic formalism for the spatial and mass distributions 
of substructures. We summarize the results in Sjl] and we 
compare our analytic predictions with numerical data in iJS] 
Finally, in ^ we discuss the results and draw final conclu- 



2 GLOBAL TIDES AND DYNAMICAL 
FRICTIONS 

2.1 Tidal Stripping 

The most dominant force that drive the dynamical evolution 
of the dark halo s ubstructures (subhalos) are the global tides 
from host halos dOkamoto fc Hab3ll999f) . The global tides 
from host halos strip the outer parts of subhalos, resulting 
in the subhalo total disruption or at least significant amount 
of subhalo mass loss. 

For the realistic treatment of the mass-loss caused by 
the global tides, it is necessary to assign the density profiles 
to both the host halos and the subhalos. We consider a sit- 
uation that a subhalo with with virial mass mvir is moving 
in a circular orbit of radius R from the center of a host halo 
with virial mass A/vir. Then we assume that the initial den- 
sity distributions of both the subhalos and the host halos 
are well described by the profiles obtai ned in N-body simu- 
lations jNavarro. Frenk fc Wliit3ll997l. hereafter NFW): 



p{R) = 



{R/R,)il + R/Rsy 



(1) 



where ps is the characteristic density that can be computed 
from the nonlinear overdensity Avir, and Ra is the halo scale 
radius. We use the top-hat radius as the virial radius such 
that 



-Rvir 



/ 3Mvir 



\ 47rAvirP 



1/3 



(2) 



where p is the mean density of the universe. The ratio of 
the scale radius Rs to the halo virial radius i?vir defines 
the halo concentrat ion parameter C and Afvir at redshift z 
jKlvpin et alJll999^ : 



C : 



124 



1 + z \ Ift-iMc 



(3) 



here we incorporate the redshift dependence found by e.g., 
iBullock et alT (1200 ill . The halo mass confined within the ra- 
dius of R can be computed from the following relation: 



M{R) = Mv 



where 



' fiC) ' 



/(X) = ln(l + X)- 



X 



1 + X' 



(4) 



(5) 



and X = R/Rs- Equations Q-JKJ are written in terms of 
the host halo properties but they also hold for the subhalo 
properties. From here on, the capital letters and the small 
letters are consistently used for the notation of host halos 
and subhalos, respectively. For example, M, R, and C repre- 
sent the mass, the radius, and the concentration parameter 
of host halos, while m, r, and c are the same quantities for 
the subhalos. 

The effect of global tide can be quantified by the tidal 
radius, rt, that is defined as the radius all the mass of a 
satellite beyond which gets lost by the global tides. For the 
realistic case that both the host halos and the subhalos are 
not point-like masses but have extended profiles as is our 
case, the tidal radius rt is given as 



n = min(rto,T-r, 



(6) 
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where rto and r^c are defined as a radius at which gravity 
equals the tidal force and a radius determined by the res- 
onances, respectively. The following two equations allow us 
to determine the values of rto and Vic'- 



fix) 
fix) 

where x = r/r.. 



an 



RaVn 



2 - 



ii+xrf{x) 



RsVn 



(7) 



(8) 



and X = R/Rs- Here Knax and Wmax repre- 
sent the maximum circular velocity of the host halo and 
the subhalo respectively. We approximate that the max- 
imum circular velo city occurs at twice the scale radius 
jKlvDin et alJll99flri . V^^^ = V{2R^) and u„,ax = v{2r,), 
where the rotation velocities are derived by 



V^{R) 



V (r) — 



R f{C) ' 

Grrivir f{x) 
r f{c) ■ 



(9) 
(10) 



Once the tidal radius is determined through the above equa- 
tions lOl- llUII . the final subhalo mass, mt, after the tidal 
stripping effect from the global tides is computed as 



rrif 



fin/rs 



(11) 



otherwise. 



fie) ' 
for rt < rvir and mf — 

2.2 Orbital Decay 

Since the subhalos reside in very dense environments within 
the host halos, they undergo dynamical friction. Although 
the most dominant force that drives the subhalo mass-loss 
is the global tides from host halos, dynamical friction also 
plays an important role in driving the subhalo dynamical 
evolution. It causes the orbital decay of the subhalos, making 
them more susceptible to strong tidal forces. 

The key quantity in describing the orbital decay caused 
by dynamical friction is the friction time scale tdi'- 

f = -^' (12) 
at tdi 

which can be estimated by the Chandraskekhar's formula: 
' d In M{R) 



tdiij. 



where 



5(C) 



In A 



,R) = o 



d\nR 



+ 1 



vlAR) 



A7vG^{lnA)m^irp{R)g{VcirciR)/V2a, 



erf (?) 



=Ce" 



» m 



2X{1 + Xf 



/(^) 



-dx. 



-(13) 

") 

(14) 
(15) 
(16) 



Here we assume that 
of the subhalos fo llow 
jBinnev fc Troma£i \l98' 



the initial circular velocities 
the Maxwellian d istributions 
iKlvpin et aP [l99i), and the 
Coulomb logarithm A has a constant value of 8 since 
this value is found to give the best result in simulations 
jTormen et alJll998h . 



Since the dynamical fricti on timescale is app roximately 
proportional to the radius R jKlvpin et al.lll999ll . the final 
decayed radius R{ of the subhalo during the time interval 
At can be estimated as 



Ri = Ri[l 



At 



tdfirriyir, Ri) 



(17) 



where Ri is the initial orbital radius of the subhalo before 
the orbital decay, and mvir is the initial virial mass of the 
subhalo before the tidal stripping effect. From here on, we 
drop the subscript "vir" in denoting the initial virial mass, 
and use the notation M and mi for representing the virial 
mass of host halos and the initial virial mass of subhalos, 
respectively. 

The model we adopt in this paper is rather simple, and 
differs from the ones used in the recent comp rehensive stud- 
of these dynamical effects. For instance, iHaveishi et alJ 



1200311 showed that the effects of tides may be underes- 
timated in the model described above, and that the im- 
pulse approximation r esults in bette r agreement with the 
numerical simulations. iBenson et alJ |2004) considered the 
improved dynamical friction model which incorporates non- 
circular motions and more complicated Coulomb logarithm 
A. Nevertheless we keep this simple model so as to be ana- 
lytically tractable. 



3 DERIVATIONS OF SUBSTRUCTURE 
DISTRIBUTIONS 

In ^ we have investigated the effects of global tides and 
dynamical frictions on the evolution of the dark matter sub- 
halos, and showed that it is possible under some simplified 
assumptions to predict the final mass and position of a sub- 
halo from the initial mass and position along with the given 
host halo mass. 

The initial spatial and mass distribution of sub- 
halos can be comput ed by modifying th e popula r 
IPress fc Schechteil Jl974l hereafre PS) approach llLeell2004ll . 
Let n{mi, Ri; M, z, Zi)dRidmi represent the number density 
of subhalos formed at redshift Zi with mass in the range of 
[mi, mi -I- drrii] located in a spherical shell of radius _Ri with 
thickness of dRi from the center of a host halo with mass 
M at redshift z. This number density can be written as a 
product of two distributions: 



n{mi,Ri;M,z,Zi) = P{Ri; M, z)n{mi\Ri; M , z, Zi), 



(18) 



where P{Ri\ M, z)dRi is the probability that the subhalo has 
an orbital radius of -Ri provided that it is included in a host 
halo of mass M at redshift z, and n(mi|iii; M, z, Zi)dmi is the 
number density of subhalos formed at Zi with initial mass in 
the range of [mi, mi -I- dmi] provided that it is located at an 
initial orbital radius of Ri from the center of a host halo of 
mass M at redshift z. 

For the probability P{Ri;M,z), we make apriori as- 
sumption that it has a form of the NFW profile, expecting 
that the subhalos follow initially the distribution of dark 
matter particles: 



P{Ri;M,z)dRi 



with 



AuRfA 



{Ri/Rs)il + Ri/R,) 



rdRi, 



(19) 
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A 



(20) 



Here the amplitude A was determined from the normaUza- 
tion constraint of 

-Rvir 

P{Ri-M,z)dRi = l. (21) 

While we normalize the distribution by equation 1211 , we ex- 
tend this distribution to _Rc , where the radius upper limit Rc 
represents the effective range of the dynamical friction force 
beyond which the force is negligible. Since the host halo has 
an extended density profile (see eq. ^), the effective range 
of dynamical friction Rc does not necessarily coincide with 
the virial radius of the host halo. Hence in order to derive 
the subhalo distribution, we have to consider those subha- 
los which were initially placed outside the virial radius of 
the host halo but eventually fell into within the host halo 
virial radius due to the orbital decay caused by the dynami- 
cal friction. We estimate the effective radius of the host halo 
dynamical friction, to find Rc ~ lOOiZa. Hereafter we use the 
value of Rc = 107?vir since i?vir ~ 107?s. As we will see in 
however, our results are insensitive to the specific choice 
of iic. 

The conditional distribution, n{mi\Ri; M, z, Zi), can be 
obtained by incorporating the spatial correlation between 
the subh alo and the host halo into the conditional PS mass 
function llYano. Naeashima fc Goudalll996l: lLe3l2004ll : 



n{mi\Ri; M, z, Zi)dmi 
where 

^_ Scizo 



2 M 



das 



drrii 



dp 



das 



Scjz) al 

^ai - at/al V 5c{z,) al 

A^(fc)dlnfc, 



f.fc(M) 

al = I A^(fc)dlnfc, 



2 / A 2 / 7 X sm kUi 

ar = 



A^(fc)- 



kRi 



dlnfc, 



e-^'/^dm;, (22) 

(23) 
(24) 
(25) 
(26) 



where A(fc) is the dimensionless power spectrum of the lin- 
ear density field, the wave number k{M) is related to the 
mass as 



k{M) = 



M 



(27) 



and 5c{z) is a threshold value of the dimensionless density 
contrast 5 at redshift z given by 5c(z) ~ 1.68/-D(z), where 
D[z) is the linear growth rate of the density field. For an 
Einstein De-Sitter universe, it reduces to 5c{z) ~ 1.68(1 + ^). 

Our use of a NFW profile form for the Lagrangian distri- 
bution, P{Ri) can be justified as follows: The original PS for- 
malism assumes that bound objects form at the local density 
peaks. Thus, in the strict P S framework, a Gaussian peak 
profile llBardeen et aljriOSfil) has to be used for P(Ri) to 
be consistent. However, sev eral numerical simulations (e.g., 
iKatz. Quinn. fc Gelblll993f l have found a poor correlation 
between the location of bound objects and the local den- 
sity peaks in Lagrangian space. In other words, in spite of 



the fact that the PS theory happens to be quite successful 
in estimating the statistically averaged number density of 
bound objects, it has been shown to fail in predicting the 
Lagrangian density profiles of bound objects, P{Ri). All one 
can say for sure is that no matter what the functional form 
of P{Ri) is, when it is mapped to the Eulerian space, it must 
be close to the NFW profile. Without knowing a correct form 
of P{Ri), a possible way to make a best approximation is 
to use the form of NFW profile itself in Lagrangian space^ . 
This guarantees that the corresponding Eulerian profile is 
NFW, and allows us to use the PS approach in Lagrangian 
space to the substructure mass function. 

The other quantity to be considered is the subhalo 
formation epoch^ distribution, dp/dz, since the orbital de- 
cay of subhal o s depe nd on the time interval At = t — U. 
iLacev fc Cold ^ll^93^ derived an a nalytic expression for 
dvldz. and'Kite^ mafc Sutd (Il996l) provided a fitting for- 
mula to dp /dz, which we adopt here (see Appendix C in 
iKitavama fc Sutdlll996l) . 

The analytic steps to reach the final subhalo spatial and 
mass distributions with taking account of tidal mass-loss and 
orbital decay caused by dynamical friction are summarized 
as follows: 

(i) We determine the final position Rf of a subhalo after 
the dynamical friction effect through equations I13ll - 117| l. 

(ii) We determine the tidal radius rt of the subhalo at the 
final position _Rf through equations lOl- llUII . 

(iii) We abandon those subhalos whose tidal radius rt is 
larger than its final orbital radius Ri, assuming that they 
will get disrupted completely by the tidal stripping effect. 

(iv) For the survived subhalos with rt < Rf, we determine 
the final mass mf through equation llll l. 

(v) Finally, we count the cumulative number of those sub- 
halos as function of mass m and radius R. 

The above procedures to evaluate the cumulative spatial and 
mass distributions of the subhalos can be summarized by the 
following analytical expression: 



N(> m,> R;M,z) = —dzi / drmdRi 
J dzi 

xn(mi, _Ri; Al, z, Zi). 
where S represents the following condition: 

mf{mi, Ri, Zi) > 



Rf{mi, Ri, z\ 

-Rv: 

-Rf (mi, Ri, Z[ 



> 
> 
> 



m, 
R, 

-Rf (mi, -Ri, Zi), 
rt{mi,Ri,Zi), 



(28) 

(29) 
(30) 
(31) 
(32) 



In other words, the integration over drmdRi is performed 
only those regions in the rrii-Ri plane that satisfy the above 



^ Since we are considering a region which has yet to virialize, it 
may have a different profile from the NFW profile; for instance, 
Sheth et al. 1 2001) considered a time-evolution of such profiles by 
modifying the scale radius. 

^ The formation epoch is defined as the epoch when a progenitor 
with a mass of fM is formed. We adopt / = 1/2 throughout this 
paper, which is a standard value used by iLacev fc Cold ilflQ.'j) 
and in many studies which is a standard choice in this context. 
We note that our result is rather insensitive to the specific choice 
of/. 
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_ 1 1 1 Mlllj 1 1 1 Mlllj 1 1 IIMLl 

: M=10i3h-iMQ : 
: m/M>10-3 : 
- z = 


_ 1 1 1 Mlllj 1 1 1 Mlllj 1 1 1 IIIU 


_ 1 1 1 Mlllj 1 1 1 Mlllj 1 1 1 IIIU 


Z;=0.2 
1 1 


Zi=1.0 

1 1 


z. = 2.0 

1 1 



10-3 10-2 10-1 lQ-3 10-2 10-1 lQ-3 lO'S 10"! 1 

m/M 

Figure 1. The region of tlie integral (eq. 1281 ') in tlie mi-Ri plane, for fixed values of zi; zi = 0.2 (left), 1.0 (middle), and 2.0 (right). 
We plot contours of the integrand (see eq. 1281 '! in the logarithmic space; (dp/dzi)n(mi, Ri; M, z, Zi)miRi. Darker contours mean larger 
values, thus darker regions contribute to the integral more. 



conditions (eqs. I29I - I32I ). To perform an integral with such a 
complicated boundary condition, it is useful to use a Monte- 
Carlo integral method. 



4 RESULT 

As a specific example, we show results of our analytic model 
in a Lambda dominated CDM model with the choice of 
the mass density parameter ilm = 0.3, the vacuum energy 
density parameter = 0.7, the spectral shape parameter 
F = 0.168, the dimensionless Hubble constant oi h = 0.7, 
and the rms fluctuation normalization erg = 0.9. We assume 
the power spectrum of cold dark matter model with primor- 
dial spectral index n j = 1, and adopt a fitting formula of 
iBardeen et all Jl986h . 

In Figure^ we show the integral region, which is deter- 
mined by equations II29|I - H32^ . in the m-^-R\ plane for fixed 
values of zi . In reality, we plot contours of the integrand (see 
eq. |28| ) in the logarithmic space. As seen, the region of the 
integral has a simple topology, and is easily understood; the 
more massive substructures tend to sink faster, thus they 
should be in the outer regions of the host halo; when Zi is 
larger the typical value of R\ also becomes larger because 
substructures experience dynamical friction more. However, 
we also find that the contribution of massive substructures 
to the integral is rather minor. This is simply because of 
the small number of massive substructures. Too small Ri is 
not allowed because final sizes of substructures should be 
smaller than the final distances from the center (eq. |32)'l. It 
is also clear from this Figure that our results are insensitive 
to the value of i?c as far as we adopt sufficiently large _Rc, 
i?c > 5/?vir, since substructures which initially lie very far 
from the host halo are not counted after all. 



TTj 1 I IllltJ 




10-6 10-6 10-4 10-3 10-a 10-1 1 
m/M 

Figure 2. Cumulative mass function of substructures in units of 
rescaled substructure mass. The redshift is z = 0. For the mass of 
the host halo, we consider M = IQl^/i-lM© (solid), IO^^Ii-'^Mq 
(dashed), and lO^/i^^M© (dotted). 



4.1 Mass Distribution 

The cumulative mass distribution A'^(> m; M, z) of subhalos 
in a given host halo of mass M at redshift z is now straight- 
forwardly obtained by setting _R = in equation 12811 : 



iV(> m;M,z) = N(> m,> Q;M,z). 



(33) 



Figure 121 plots A'^(> m;M,z) at redshift 2; = as a 
function of rescaled subhalo mass for the three different cases 
of host halo mass M. Note that all the three curves in Figure 
|21are very similar to one another, indicating that the subhalo 
mass distribution is almost independent of the host halo 
mass. In addition, these curves are close to a power law with 
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Figure 3. Normalized cumulative spatial distribution of substructures in units of rescaled radius. Mass of the host halo is set to 
M = lO^^/i~lM0 (left), lO^^h~^M0 (middle), and lO^^/i^^M© (right). We only consider substructures with mass larger than 10~^M 
(solid), 10~*M (dashed), and 10~^M (dotted). The spatial distribution of smooth dark matter component, which can be calculated from 
equation 0, is also shown by dash-dotted lines. 



the slope of ~ —0.9 in low mass range of m/M ~ 10 ^ and 
~ — 1.0 for high mass ranee of m/M ~ 10-2. These findings 
from our analytic model are all consistent wit h recent high- 
resolution numerical simul ation results (e.g., iGhigna et alJ 
I2OOOI: iDe Lucia et aLlliool . It is worth noticing, however, 
that there is some tendency that the larger the host halos 
are the slightly more massive subhalos they are assigned to, 
which is likely because larger host halos form relatively late 
so that they have relatively little time for their subhalos to 
undergo orbital decay and lose mass. 



4.2 Spatial Distribution 

Equation 12811 can be interpreted as the cumulative spatial 
distribution of subhalos for a given mass range. Figure |3 
plots the normalized cumulative spatial distribution of sub- 
halos as a function of rescaled radius, and the spatial distri- 
butions of the smooth component dark matter for compar- 
ison as well. The normalized spatial distribution is defined 
as N{> m, > R; M, z)/N{> m; M, z). One can see that the 
spatial distributions of the dark subhalos are anti-biased rel- 
ative to that of the dark matter particle components. This 
phenomenon has been observed in several numerical simu- 
lations (e.g.. lGhigna'et aL .2000: De Lucia et al. 20ol) . We 
also note that more massive subhalos have stronger tendency 
of anti-bias, which is likely because this anti-bias is caused 
by tidal stripping and dynamical friction, so that massive 
subhalos are more affected by those dynamical processes. 



4.3 Velocity Distribution 

The circular velocity distribution of subhalos is sometimes 
more useful in practice bec ause it is more readil y to be com- 
pared with observations llCole fc Kaiseil |l989l JShunas^a 
I1993I: iGonzalez etHI I2OO0I: ISheth et alJ l200a iDesai et al " 
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Figure 4. Cumulative circular velocity function of substructures 
in units of rescaled substructure velocity. For the mass of the 
host halo, we consider M = IO^^H-'^Mq (solid), IO^^H-^Mq 
(dashed), and 10^^ H'^Mq (dotted). 



120041) . Our analytic model also enables us to compute the 
circular velocity distribution. Without the tidal stripping ef- 
fect, the subhalo would retain the original NFW density pro- 
file, and thus the subhalo circular velocity would be given by 
«max ~ v{2rs). However, because of the tidal stripping effect, 
the subhalos eventually end up with having truncated den- 
sity profile at truncation radius ~ 2rs (Havashi et al. 200,^. 
Thus, in our realistic model, the subhalo circular velocity is 
given as 



'^circ — 



v(2rs; mi, Zi) 
v(rt;mi,Zi) 



in > 2rs), 
(r-t < 2rs). 



(34) 
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Thus, the subhalo circular velocity depends not only on the 
subhalo mass but also on the subhalo position and formation 
epoch. Using this relation (eq. |84| 'l along with the subhalo 
spatial and mass distribution (eq. 1281 ). we can also evaluate 
the subhalo velocity distribution. 

Figure 13 plots the cumulative circular velocity distribu- 
tion of subhalos versus the subhalo velocity rescaled by Knax. 
Note that the subhalo velocity distribution also shows very 
weak dependence on the host halo mass, power-law shaped 
with the slope of ~ — 3 around Ucirc/Knax ~ 10~^, which 
is all consistent with numerical detections dOhiena et alJ 
|200(J) . And it also has the same tendency as the spatial 
distribution that the larger a host halo is the more subhalos 
with high velocity it has. 



5 COMPARISON WITH NUMERICAL 
SIMULATIONS 

In this section, we compare our analytic predictions with 
nu merical data from high -resolution simulations performed 
bv lDe Lucia et all ||2004) . In addition to our own analytic 
model, we also compare the following two models with nu- 
merical simulations to demonstrate the importance of dy- 
namical processes that we have included: (i) the model with 
tidal stripping but without dynamical friction, (ii) the model 
without tidal stripping and dynamical friction. The models 
(i) an d (ii) basically correspon d to those considered in iLed 
l)2no4) andlFuiita et alJ J2Qq3)) respectively. 

First in Figure |3 we plot our analytic predictions of 
the subhalo mass distribut ion with the numeric al simula- 
tion results (see Fig. le of iDe Lucia et al]|2004) . We find 
that our model agrees quite well with the numerical simula- 
tions, while the models without dynamical friction or tidal 
stripping effect tend to overpredict the number of massive 
substructures. Those models are inconsistent with the nu- 
merical simulations even if we allow them to change their 
normalization. The success of our model indicates that the 
effects of tidal stripping and dynamical friction are essen- 
tial to understand the evolution of subhalos and understand 
theh distributions. For M = IO^^/i'^Mq and lO"/i"^M0, it 
seems that there is a small difference at m/M > 0.02. How- 
ever, this difference is not significant given the small num- 
bers of host halos (5 and 6, respectively) used for studying 
the mass distribution in the numerical simulations. 

In Figure |S| we compare our analytic predictions on the 
subhalo spatial distribution with numerical data (see Figure 
7 of Dc Lucia ct al. 2004), and show that the two results 
agree with each other quite well. It is clear from this Figure 
that our model successfully reproduces the feature found in 
the numerical simulations that more massive substructures 
are preferentially located in the outer part of their host halo. 
This is not the case with the model without dynamical fric- 
tion. The model without tidal stripping and dynamical fric- 
tion also shows this feature, because we have assumed also 
in this model that those subhalos whose radius rvir is larger 
than its final orbital radius R{ will get disrupted completely. 
However, the distributions themselves are largely inconsis- 
tent with the numerical distributions, and rather similar to 
the distribution of dark matter particles as for smaller sub- 
structures. It implies that the dynamical frictions is the main 



reason for the anti-bias of the subhalo spatial distribution 
relative to that of the dark matter particles. 



6 DISCUSSIONS AND CONCLUSIONS 

We have constructed an analytic model for mass and spa- 
tial distributions of dark subhalos by taking account of two 
dominant dynamical processes that drive dominantly sub- 
halo evolution: one is tidal stripping and consequent mass- 
loss caused by the host halo global tides which remove the 
outer, and the other is the orbital decay caused by dynami- 
cal friction which drives massive subhalos to the inner part 
of the host halo. We also incorporate the formation epoch 
variation of the host halo, and the orbital decay of satellite 
halos outside the host halo virial radius. 

We have found that our model predicts nearly power- 
law mass distribution with weak dependence of host halo 
mass. The slope of the cumulative mass function turned out 
to be about —0.9 for low mass substructures m/M ~ 10~^, 
and about —1.0 for larger mass substructures m/M ~ 10~^. 
Next we have predicted the spatial distributions for given 
mass ranges of substructures. We have found that the spatial 
distributions are basically anti-biased relative to the smooth 
component, i.e., the substructures are preferentially located 
in the outer part of the host halo. The extent of the anti- 
bias is larger for more massive subhalos. Using our model, 
we have also predicted the velocity distributions of subhalos. 
These findings are all consistent with what has been detected 
by recent high-resolution numerical simulation. So, we have 
provided physical explanations to those characteristic find- 
ings of numerical simulations on the subhalo distributions. 

There is some discussion about the accuracy of the 
iBullock et alJ i200ll) choice of th e evolution of the concen- 
tration parameter. For instance, IZhao et alJ (|200.?l ) found 
that all high-redshift massive halos have a similar median 
concentration, C ~ 3.5. To see how this affects our results, 
we repeat the computation by assuming the concentration 
parameter of equation Q if C ^ 3.5, and replace it with 
C — 3.5 if equation Q yields C < 3.5. We find that the 
substructure mass functions are not changed with this re- 
placement; at 2: < 1 the change is unnoticeable, and even at 
2: = 2 it changes the substructure mass function by < 10%. 

To test the validity of our model, we have compared 
our analytic model with the results of th e high-resolution 
numerical simulations by iDe Lucia et alJ li2004il . We have 
showed that mass and spatial distributions of substructures 
in our model excellently agrees with those in the numerical 
simulations. For comparison, we have considered the mod- 
els without tidal stripping or dynamical friction and found 
that they cannot reproduce the distributions in the numeri- 
cal simulations. This indicates that the effects of both tidal 
stripping and dynamical friction are very essential for con- 
structing a realistic model of subhalo distributions. 

Yet, it should be noted that we have adopted several 
simplified assumptions. For instance, we have ignored the 
effect of encounters between substructures, although it may 
lead to mass loss comparable to that caused by global tidal 
stripping ( Tormen et al. 1 998) . We have also assumed spher- 
ical halos and circular orbits of subhalos, which are in- 
accurate because halos in the CDM model are shown to 
be rather triaxial jjing fc Sutdl2002h and the subhalo or- 
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Figure 5. Comparis on of the analytic mass distributions {thick lines) witli numerical ones {thin lines and triangles) computed by 
iDe Lucia et alJ (20^. We consider the host halo mass of M = W^^H-'^Mq {upper left), lO^^/i'^Mo {upper right), lO^^h'^M© {lower 
left), and lQ^^h~^ Mq {lower right). As for the analytic mass distributions, we consider not only the fiducial mass distributions in which 
both dynamical friction and tidal stripping are included {solid), but also the mass distributions without dynamical fricti on (dashed) and 
witho ut dynamical friction and tidal stripping (dotted). Here we use N{m) + 1 rather than N(> m) because in Figure le of lDe Lucia et all 
i200'J) the host halo itself was included in the cumulative mass function. 




Figure 6. Comparis on of the analytic spatial distributions (thick lines) with numerical ones (thin lines and symbols) computed by 
iDe Lucia et alJ i2004h . Comparison is done for both massive substructures (m > 10~^M, denoted by solid lines and/or crosses) and less 
massive substructures (10~'^M > m > 4 x 10~*M, denoted by dashed lines and/or triangles). From left to right panels, we compare the 
fiducial analytic model, the model without dynamical friction, and the model without dynamical friction and tidal stripping. 
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bits are shown to be quite eccentric jTormen et al.lll998t 
iHavashi et al]l2003ri . In spite of these drawbacks of our an- 
alytic model, we believe that our model for the subhalo 
distributions is the most realistic one that has ever been 
suggested, and that it can be applied to wide fields such 
as gravitational lensing and halo model of large-scale struc- 
ture. Our analytic model is quite useful also in calculating 
the mass and spatial distributions of substructures for dif- 
ferent initial conditions such as running primordial power 
spectrum suggested by the combined analysis of cosmic mi- 
crowave background and large-scale structure observations 
llSoergel et alJl200.1 ). 
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